library(tidyverse)
library(haven)
library(stargazer)
library(broom)
library(scales)
library(margins)
library(kableExtra)

#### Study 1 ####
#### Load Lucid sample ####
df <- read.csv("climate_lucid.csv")
## Local news readership
df <- df %>% mutate(local.read = case_when(Q10.8 == "Never" ~ 1,
                                           Q10.8 == "Occasionally" ~ 2,
                                           Q10.8 == "Several times a week" ~ 3,
                                           Q10.8 == "Everyday" ~ 4))
df <- df %>% mutate(national.read = case_when(Q10.9 == "Never" ~ 1,
                                              Q10.9 == "Occasionally" ~ 2,
                                              Q10.9 == "Several times a week" ~ 3,
                                              Q10.9 == "Everyday" ~ 4))



## manipulation checks
df <- df %>% mutate(pass1 = case_when(Q2.2=="The Advocate" ~ 1,
                                      Q2.2=="New York Times" ~ 0,
                                      Q3.2=="The Advocate" ~ 1,
                                      Q3.2=="New York Times" ~ 0,
                                      Q4.2=="The Advocate" ~ 1,
                                      Q4.2=="New York Times" ~ 0,
                                      Q5.2=="The Advocate" ~ 0,
                                      Q5.2=="New York Times" ~ 1,
                                      Q6.2=="The Advocate" ~ 0,
                                      Q6.2=="New York Times" ~ 1,
                                      Q7.2=="The Advocate" ~ 0,
                                      Q7.2=="New York Times" ~ 1
))

## Local news readership
df <- df %>% mutate(local.read = case_when(Q10.8 == "Never" ~ 1,
                                           Q10.8 == "Occasionally" ~ 2,
                                           Q10.8 == "Several times a week" ~ 3,
                                           Q10.8 == "Everyday" ~ 4))
df <- df %>% mutate(national.read = case_when(Q10.9 == "Never" ~ 1,
                                              Q10.9 == "Occasionally" ~ 2,
                                              Q10.9 == "Several times a week" ~ 3,
                                              Q10.9 == "Everyday" ~ 4))

df <- df %>% mutate(female = case_when(Q10.3=="2.\tFemale" ~ 1,
                                       Q10.3=="1.\tMale" ~ 0,
                                       Q10.3=="3. Other" ~ 0),
                    age = Q10.2_1,
                    income = case_when(Q10.4=="Less than $10,000" ~ 1,
                                       Q10.4=="$10,000 - $19,999" ~ 2,
                                       Q10.4=="$20,000 - $29,999" ~ 3,
                                       Q10.4=="$30,000 - $39,999" ~ 4,
                                       Q10.4=="$40,000 - $49,999" ~ 5,
                                       Q10.4=="$50,000 - $59,999" ~ 6,
                                       Q10.4=="$60,000 - $69,999" ~ 7,
                                       Q10.4=="$70,000 - $79,999" ~ 8,
                                       Q10.4=="$80,000 - $99,999" ~ 9,
                                       Q10.4=="$100,000 - $119,999" ~ 10,
                                       Q10.4=="$120,000 - $149,999" ~ 11,
                                       Q10.4=="$150,000 or more " ~ 12),
                    white = case_when(Q10.5=="White" ~ 1, 
                                      Q10.5=="Black" ~ 0,
                                      Q10.5=="Hispanic" ~ 0,
                                      Q10.5=="Asian" ~ 0,
                                      Q10.5=="Middle eastern" ~ 0,
                                      Q10.5=="Mixed" ~ 0,
                                      Q10.5=="Native American" ~ 0,
                                      Q10.5=="Other" ~ 0))
#### Descriptive Stats of Covariates (Table S6) ####
df$age <- as.numeric(df$age)
summary(df[, c("female", "age", "income", "white", "local.read", "national.read")])
apply(df[, c("female", "age", "income", "white", "local.read", "national.read")], MARGIN=2, FUN=function(x) sd(x, na.rm=T))

#### Descriptive Stats of Outcome Variables by Party ID #### (Table S7)
df %>% group_by(pid3) %>% dplyr::select(c("Q8.2_1", "Q8.4_1", "serious", "action")) %>%
  summarise(across(where(is.numeric), mean, na.rm=TRUE)) %>% kbl(format="latex")

df %>% group_by(pid3) %>% dplyr::select(c("Q8.2_1", "Q8.4_1", "serious", "action")) %>%
  summarise(across(where(is.numeric),sd, na.rm=TRUE)) %>% kbl(format="latex")

#### Balance Tests (Table S10) ####
df$pid3 <- relevel(as.factor(df$pid3), ref="Ind")

m1<-(glm(local ~ as.factor(pid3) + female + as.numeric(age) + income + white + local.read + national.read,data=df,family=binomial(link = "logit")))
stargazer(m1)

#### Analysis Controlling for Covariates (Table S12) ####
mod1 <- lm(Q8.2_1 ~ local+ female + income +  white +  local.read + national.read, subset(df, pid3=="Rep"))
mod2 <- lm(Q8.4_1 ~ local+ female + income + white + local.read + national.read, subset(df, pid3=="Rep"))
mod3 <- lm(serious ~ local+ female + income + white + local.read + national.read, subset(df, pid3=="Rep"))
mod4 <- lm(action ~ local+ female + income + white + local.read + national.read, subset(df, pid3=="Rep"))
stargazer(mod1, mod2, mod3, mod4,
          type = "latex", header=FALSE, title = "Local Treatment Effects",
          dep.var.labels = c( "Accurate","Relevant","Seriousness", "Action"), omit.stat=c("f", "ser"))




#### Subgroup by Coastal Residence (Table S15)####
library(zipcodeR)
zipcode.dist <- table(df$Q10.6[grep("[[:digit:]]", df$Q10.6)])
zipcode.dist <- as.data.frame(zipcode.dist)
names(zipcode.dist)[1] <- "zipcode"
zipcode.dist[,1] <- as.character(zipcode.dist[,1])
zipcodes <- sort(df$Q10.6[grep("[[:digit:]]", df$Q10.6)])
zipcodes <- geocode_zip(zipcodes)
zipcodes <- zipcodes[!is.na(zipcodes$lat),]
zipcodes <- merge(zipcodes, zipcode.dist, by="zipcode")

df$zipcode <- as.character(NA)
df[grep("[[:digit:]]", df$Q10.6), "zipcode"]  <- df[grep("[[:digit:]]", df$Q10.6), "Q10.6"]

zipcodes$county <- sapply(zipcodes[,1], FUN = function(x) reverse_zipcode(x)$county)
df <- merge(df, zipcodes, by="zipcode", all.x=TRUE)
coast <- c("Cameron Parish", "Iberia Parish", "Jefferson Parish", "Lafourche Parish", "Orleans Parish", "Plaquemines Parish", "St. Bernard Parish", "St. Mary Parish", "St. Tammany Parish", "Terrebonne Parish", "Vermilion Parish")
df <- df %>% mutate(coast = ifelse(county %in% coast, 1, 0)) 
mod1 <- lm(Q8.2_1 ~ local * coast, subset(df, pid3=="Rep"))
mod2 <- lm(Q8.4_1 ~ local * coast, subset(df, pid3=="Rep"))
mod3 <- lm(serious ~ local * coast, subset(df, pid3=="Rep"))
mod4 <- lm(action ~ local * coast, subset(df, pid3=="Rep"))
stargazer(mod1, mod2, mod3, mod4, 
          type = "latex", header=FALSE, title = "Perceptions of Article",
          dep.var.labels = c( "Accurate", "Relevant","Serious", "Action"), omit.stat=c("f", "ser"))

#### Study 2 ####
df <- read.csv("climate_qualtrics.csv")

#### Descriptive Statistics ####
df <- df %>% rename(age = Q16.2_1)
df$age <- as.numeric(df$age)
df <- df %>% mutate(female = case_when(Q16.3=="Female" ~ 1,
                                       Q16.3=="Male" ~ 0,
                                       Q16.3=="Other" ~ 0),
                    income = case_when(Q16.4=="Less than $10,000" ~ 1,
                                       Q16.4=="$10,000 - $19,999" ~ 2,
                                       Q16.4=="$20,000 - $29,999" ~ 3,
                                       Q16.4=="$30,000 - $39,999" ~ 4,
                                       Q16.4=="$40,000 - $49,999" ~ 5,
                                       Q16.4=="$50,000 - $59,999" ~ 6,
                                       Q16.4=="$60,000 - $69,999" ~ 7,
                                       Q16.4=="$70,000 - $79,999" ~ 8,
                                       Q16.4=="$80,000 - $99,999" ~ 9,
                                       Q16.4=="$100,000 - $119,999" ~ 10,
                                       Q16.4=="$120,000 - $149,999" ~ 11,
                                       Q16.4=="$150,000 or more " ~ 12),
                    white = case_when(Q16.5=="White" ~ 1, 
                                      Q16.5=="Black" ~ 0,
                                      Q16.5=="Hispanic" ~ 0,
                                      Q16.5=="Asian" ~ 0,
                                      Q16.5=="Middle eastern" ~ 0,
                                      Q16.5=="Mixed" ~ 0,
                                      Q16.5=="Native American" ~ 0,
                                      Q16.5=="Other" ~ 0),
                    college = case_when(Q16.6 == "4-year College Degree" ~ 1,
                                        Q16.6 == "Graduate or professional degree" ~ 1,
                                        Q16.6 == "Some College" ~ 1,
                                        Q16.6 == "High School" ~ 0,
                                        Q16.6 == "Less than High School" ~ 0),
                    local.read =  case_when(Q16.12 == "Never" ~ 1,
                                            Q16.12 == "Occasionally" ~ 2,
                                            Q16.12 == "Several times a week" ~ 3,
                                            Q16.12 == "Everyday" ~ 4),
                    nyt.read = case_when(Q16.13_1 == "Never" ~ 1,
                                              Q16.13_1 == "Occasionally" ~ 2,
                                              Q16.13_1 == "Several times a week" ~ 3,
                                              Q16.13_1 == "Everyday" ~ 4),
                    usa.read = case_when(Q16.13_3 == "Never" ~ 1,
                                         Q16.13_3 == "Occasionally" ~ 2,
                                         Q16.13_3 == "Several times a week" ~ 3,
                                         Q16.13_3 == "Everyday" ~ 4))
#### Table S8 ####
summary(df[, c("female", "age", "income", "white", "local.read", "nyt.read", "usa.read")])
apply(df[, c("female", "age", "income", "white", "local.read", "nyt.read", "usa.read")], MARGIN=2, FUN=function(x) sd(x, na.rm=T))


#### Descriptive Stats of Outcome Variables by Party ID #### (Table S9)
df <- df %>% mutate(repub = ifelse(Q3.1 == "Republican" | Q3.4 == "Republican Party", 1, 0))
df %>% group_by(repub) %>% dplyr::select(c("credible", "relevant", "happening", "serious", "harm.when", "action")) %>% 
  summarise(across(where(is.numeric), mean, na.rm=TRUE)) %>% kbl(format="latex")

df %>% group_by(repub) %>% dplyr::select(c("credible", "relevant", "happening", "serious", "harm.when", "action")) %>% 
  summarise(across(where(is.numeric), sd, na.rm=TRUE)) %>% kbl(format="latex")

#### Balance Test (Table S11) ####
m1<-(glm(as.factor(FL_8_DO == "LocalNews,Treatment") ~  female + as.numeric(age) + income + white + college,data=df,
         family=binomial(link = "logit")))
stargazer(m1)

#### Analysis Controlling for Covariates (Table S13)####
mod1 <- lm(credible ~ as.factor(FL_8_DO) + female + age +  white + college+ local.read , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod2 <- lm(relevant ~ as.factor(FL_8_DO) + female +  age + white + college+ local.read, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod3 <- lm(happening~ as.factor(FL_8_DO) + female + age +  white + college+ local.read , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod4 <- lm(serious ~ as.factor(FL_8_DO)+ female +  white + college+ local.read , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod5 <- lm(harm.when ~ as.factor(FL_8_DO)+ female +  age + white + college+ local.read, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican"| Q3.4 == "Republican Party"))
mod6 <- lm(action ~ as.factor(FL_8_DO) + female +  age + white +  college + local.read, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
stargazer(mod1, mod2, mod3, mod4, mod5, mod6,
          type = "latex", header=FALSE, title = "Local Treatment Effects",
          dep.var.labels = c( "Credible" ,"Relevant", "Happening", "Serious", "Harm", "Action"), omit.stat=c("f", "ser"))

## Strong Republicans (Table S14)
mod1 <- lm(credible ~ as.factor(FL_8_DO) + female + age +  white + college+ local.read , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" & Q3.2 == "Strong"))
mod2 <- lm(relevant ~ as.factor(FL_8_DO) + female +  age + white + college+ local.read, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" & Q3.2 == "Strong"))
mod3 <- lm(happening~ as.factor(FL_8_DO) + female + age +  white + college+ local.read , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" & Q3.2 == "Strong"))
mod4 <- lm(serious ~ as.factor(FL_8_DO)+ female +  white + college+ local.read , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" & Q3.2 == "Strong"))
mod5 <- lm(harm.when ~ as.factor(FL_8_DO)+ female +  age + white + college+ local.read, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican"& Q3.2 == "Strong"))
mod6 <- lm(action ~ as.factor(FL_8_DO) + female +  age + white +  college + local.read, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" & Q3.2 == "Strong"))

stargazer(mod1, mod2, mod3, mod4, mod5, mod6,
          type = "latex", header=FALSE, title = "Local Treatment Effects",
          dep.var.labels = c( "Credible" ,"Relevant", "Happening", "Serious", "Harm", "Action"), omit.stat=c("f", "ser"))

#### Coastal Residence (Table S16) ####
library(zipcodeR)
zipcode.dist <- table(df$Q16.7[grep("[[:digit:]]", df$Q16.7)])
zipcode.dist <- as.data.frame(zipcode.dist)
names(zipcode.dist)[1] <- "zipcode"
zipcode.dist[,1] <- as.character(zipcode.dist[,1])

zipcodes <- sort(df$Q16.7[grep("[[:digit:]]", df$Q16.7)])
zipcodes <- geocode_zip(zipcodes)
zipcodes <- zipcodes[!is.na(zipcodes$lat),]
zipcodes <- merge(zipcodes, zipcode.dist, by="zipcode")

zipcodes$county <- sapply(zipcodes[,1], FUN = function(x) reverse_zipcode(x)$county)


df$zipcode <- as.character(NA)
df[grep("[[:digit:]]", df$Q16.7), "zipcode"]  <- df[grep("[[:digit:]]", df$Q16.7), "Q16.7"]
df <- merge(df, zipcodes, by="zipcode", all.x=TRUE)

coast <- c("Cameron Parish", "Iberia Parish", "Jefferson Parish", "Lafourche Parish", "Orleans Parish", "Plaquemines Parish", "St. Bernard Parish", "St. Mary Parish", "St. Tammany Parish", "Terrebonne Parish", "Vermilion Parish")
df <- df %>% mutate(coast = ifelse(county %in% coast, 1, 0)) 


mod1 <- lm(credible ~ local, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod2 <- lm(relevant ~ local, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod3 <- lm(happening~ local , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))

mod4 <- lm(serious ~ local, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod5 <- lm(harm.when ~ local, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican"| Q3.4 == "Republican Party"))
mod6 <- lm(action ~ local, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))


mod1.1 <- lm(credible ~ local * coast, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod2.1 <- lm(relevant ~ local * coast, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod3.1 <- lm(happening~ local * coast , data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))

mod4.1 <- lm(serious ~ local * coast, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))
mod5.1 <- lm(harm.when ~ local * coast, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican"| Q3.4 == "Republican Party"))
mod6.1 <- lm(action ~ local * coast, data = df %>% filter(treatment==1) %>%
             filter(Q3.1 == "Republican" | Q3.4 == "Republican Party"))

stargazer(mod1.1, mod2.1, mod3.1, mod4.1, mod5.1, mod6.1,
          type = "latex", header=FALSE, title = "",omit.stat=c("f", "ser"))
